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Abstract: This paper deals with the direction-of-arrival (DOA) estimation of local scattered 
code-division multiple access (CDMA) signals based on a particle swarm optimization 
(PSO) search. For conventional spectral searching estimators with local scattering, the 
searching complexity and estimating accuracy strictly depend on the number of search grids 
used during the search. In order to obtain high-resolution and accurate DOA estimation, a 
smaller grid size is needed. This is time consuming and it is unclear how to determine the 
required number of search grids. In this paper, a modified PSO is presented to reduce the 
required search grids for the conventional spectral searching estimator with the effects of 
local scattering. Finally, several computer simulations are provided for illustration 
and comparison. 

Keywords: particle swarm optimization; direction-of-arrival estimation; CDMA; local 
scattering 



1. Introduction 

Adaptive array techniques have been developed for enhancing the performance of code-division 
multiple access (CDMA) systems. In the CDMA system, each user employs a unique pseudo-noise 
(PN) codeword to identify their data stream. Each user's transmission interferes with all the other users 
and causes multiple access interference (MAI). The use of antenna arrays as a tool for improving 
coverage, reducing interference, and increasing capacity in CDMA systems has attracted significant 
interest [1]. Multiple propagation is common in cellular systems, which may result from local scatters 
in the vicinity of the sources. In macrocellular environments, high base station antennas will receive 
these locally scattered signals from the mobile terminal, which are coherent and confined to a narrow 
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angular region [2]. In addition, the observation period is assumed to be short in comparison to the 
coherence time of the channel so that the channel may be modeled as time-invariant. A mobile terminal 
with local scattering can be well approximated by a single point source, as seen from the base station. 
The model can express the mobile and coherent interference from the surrounding environment under 
non line-of-sight conditions as another arithmetical model. In cellular radio systems, antenna array 
processing is considered to be a useful method to solve problems such as multipath and co-channel 
interference and increase capacity [1]. There have been some studies appraising the impact of a local 
scattering channel model on CDMA system, which propose uniform linear array (ULA) geometry [3]. 
They use a first-order Taylor expansion to linearize the nonlinear conventional array manifold 
produced by scatters and develop the generalized array manifold (GAM) model which can obtain 
better nominal DOA estimation for mobiles. Here, the multiple signal classification (MUSIC) 
searching function derived by the steering vector of the GAM model is termed as GMUSIC. However, 
in all these studies, estimation of direction of arrival (DOA) of the desired signals is required. Array 
outputs aligned with code-matched filter can make the DOA estimation of multiple sources 
equivalent to that for single source localization problem in a noisy environment. With the advantage 
of code-matched filter inherent in the CDMA system, it has been proved that the multiple 
signal classification estimator [4] can obtain an unbiased DOA estimation with low mean-square-error 
(MSE) [5]. It also contributes to solving the limitation that the number of array elements must be more 
than the number of impinging sources, but for a conventional spectral searching DOA estimator such 
as MUSIC, its searching complexity and estimating accuracy strictly depend on the number of search 
grids used during the search. This is time consuming and it is unclear how to determine the required 
number of search grids. Thus, in this paper, we will focus on the GAM model of ULA to employ the 
DOA estimation in multipath reflection scenario from scatters. 

Particle swarm optimization (PSO) is a population-based stochastic optimization paradigm, in 
which each agent, named particle, of the population, named swarm, is thought of as a collision-proof 
bird and used to represent a potential solution [6]. Like a genetic algorithm (GA), PSO starts by 
initializing a population of random solutions and searches for optima by updating generations. But, 
PSO does not use any evolution operators. In PSO, particles fly through the problem space by 
following their own experience and the best experience attained by the swarm as a whole. In contrast 
to analytical or general heuristic methods, PSO is computationally efficient and has great capability of 
escaping local optima [7,8]. In addition, a key characteristic of PSO is that the algorithm itself is 
highly robust yet remarkably simple to implement, while processing similar capabilities as other 
evolutionary algorithms such as GA [9]. A maximum likelihood (ML) criterion with hard-constraint 
PSO based solution applied to DOA estimation is proposed in [10], that explores the potentially 
superior performance at less computational costs. The same PSO algorithm based on the ML 
methodology is also derived in [11], where the cost function is an extension of the ML criteria that 
were originally developed for angle estimation with some of the sensors that are perfectly calibrated, 
while others are uncalibrated. Recently an adaptive PSO algorithm has been successfully applied in 
ML optimization solutions [12]. In another report a combined fuzzy adaptive PSO and differential 
evolution are also used to solve economic dispatch problems [13]. In this paper, a modified adaptive 
inertia weight PSO (APSO) is proposed to rationally balance the global exploration and local 
exploitation abilities of the particle to the GMUSIC criteria function [3] for CDMA signals. The 
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resulting estimator is called APSO-GMUSIC, where a simple and effective measure, individual search 
ability, is defined to indicate whether each particle lacks global exploration ability in each dimension. 
By employing APSO to treat our optimization problem of the desired signal direction angle 6 . The 
DOA, obtained by the APSO-GMUSIC estimator, converges to the optimal or near optimal solution. 
Simulation results show that the proposed APSO-GMUSIC estimator is very suitable to treat the DOA 
estimation for CDMA signals in a local scattered scenario. 

2. System Model 

2.1. Array Data Model 

Consider a DOA estimation scenario in a CDMA system with P users. The data bit b p {t) e {— 1, 1} 
of the pth user is spread by a pseudo-noise (PN) signature waveform c p (t) for p = 1, 2, P. The 
signature waveform of the pth user is composed of a spreading sequence of L chips, i.e., 
c p = S P iPt l {t~(l~ l)^ c ] > where g p l e {-II VZ, +11 VZ} is the /th spreading chip for the pth 
user and prdt) is the spreading pulse of duration T c . Thus, the transmitted signal of the pth user during 
a bit interval Tb can be represented by: 



where e p is the signal amplitude. Let the processing gain of the spreading code be L = TblT c . To 
simplify analysis, the PN code during a bit interval at a sample rate \IT C can be expressed as an L x 1 
vector of discrete sequence c p = [c p i, c P 2,..., c p l] T where the superscript Tis the transpose operator. 

Employing the code -matched filter and DOA estimation algorithm, we propose a receiving base 
station of CDMA communication system with a uniform linear array (ULA). The antenna array 
consists of M identical omidirectional elements and receives line-of-sight signals from P users that 
arrive at the array from different bearing angles d\, 62, . . ., 0 P with respect to the broadside of array in a 
cell/sector. Assume that the array element space is d. The direction vector associated with the pth user 
is given by: 



where a m (0 P ) = exp[-j7td(m-\)sm(0 p )/ ft] denotes the response of the mth sensor array to a signal with 
unit amplitude arriving from the direction angle 0 P , where /? is the propagation speed of plane wave. 
Thus, the received baseband signal across the array can be written in vector form: 



s p (t) = £ p b p (t)Cp(t) 



(1) 



a(0p) = [a\(6 p ), a 2 (0 p ), aM(0 p )] T 



(2) 



p 




(3) 



where s p (t)= ^£ p b p (i)c p (t-iT h -T p ), x p is the time delay, and the observation noise vector n(t) 



i=- 



represents the spatially and temporally zero mean white Gaussian noise vector. 
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2.2. MUSIC Algorithms for Local Scattering Channels 

To show how the spatial signatures depend on the spatial distribution of the multipath propagation, 
we assume that the time dispersion produced by the multipath propagation is small compared to the 
reciprocal of the signal's bandwidth, the time delay may be approximated as a phase shift 
s p (t-T ph ) = e- J nhT * h s p (f), where/, is the carrier frequency, and x p h is the time delay associated with 

the /zth scattered signal from the pth source arriving at the array. Let s p (t) be the signal scattered by the 
pth user. Due to the multipath propagation near the pth user, its contribution to the array is modeled as 
a superposition of N p scattered rays. Then, the data received at the array is given by M x 1 vector form: 

p N k 

x(o = £ X P P h s P ^ 2 * fcTph *P> &p + ~ A Ph ) + n (o (4) 

p=\ 4=1 

where the reflection coefficient f} p h is the complex amplitude of the /zth scattered signal from the pth 
user, a. p h is the response vector of the /zth scattered signal from the pth user, and j = 4-1 . A ph is 

spread angle of the pth user due to local scattering. It is assumed that the power is normalized so that all 
rays have equal power \/N p . Then fi ph - e' Xph j ^N~ p , where Xph is i.i.d. and uniform over [0, 2;r] [3]. 

Then, the received data vector \(t) can be expressed as: 

x(0 = Us(0 + n(0 (5) 
where U = [ui, 112, u p ] and s(t) = [s\(t), sjit), s p (t)] T . The pth column of U, which is denoted as 

h=l 7 P h a P h and 

Y P h ~ fi P h e j2KfcT " h . In order to pick out the pth user's signal, a code-matched filter containing the 

specified PN code vector c p is applied to x(/). For simplicity, the bit sampling index is dropped to make 
the equations more readable in the following description. Therefore, the pth user's despread and 
sampled array vector signal y p can be represented as: 

p 

y p =xc p =Ls p u(0 p )+ X s q u(0 q ) + n p (6) 

?=i, q*p 

where s p = e p b p and s q = e p b p K qp . Hence, K =E{\s p | 2 } and K =E{\s q | 2 } denote the average 
power of the desired pth user and each interferer user, respectively. The notation E{»} is used to 
denote the expectation operator. K qp is the inner product of two different PN code vectors c q T c p , and the 
M x 1 vector n p with zero mean and variance n n is the code-matched filter output of the noise vector n. 
It was shown in Appendix A of [3] that K qp approaches a real Gaussian random variable with zero mean 
and unit variance for the normalized PN codes. Consequently, the signals of the undesired users 
through the code-matched filter, i.e., MAI, can be viewed as a noise vector n M Ai P with zero mean and 
variance tzmai p and n MA1 = x ^ p s q u(0 q ) . The correlation matrix of Equation (6) in the observation 

interval is given by: 

=E{y p y H p } = L\u(0 p )u H (d p ) + x np I M (7) 
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where the subscript H denotes conjugate transpose, n np = kmai p + n n , and Im is an MxM 
identity matrix. 

For high angular resolution subspace methods, the MUSIC algorithm is a kind of DOA estimation 
technique based on eigenvalue decomposition (EVD), which is also called the noise subspace-based 
method. The eigendecomposition of Equation (7) can be expressed as: 

M 

R P = I^e^e* = V/! + E > A „» E J' (8) 

m=\ 

where rj pl > ?] p2 = rj p3 =■■■ = jj pM = n np are the eigenvalues of R p and e pm denotes the eigenvector 
associated with J] pm for m = 1, 2, ...,M. Moreover, e p \ and E pn - [e p2 , e p3 , ••• , e pM ] are orthogonal 
and span the signal and noise subspace corresponding to K p , respectively. A pn =ft„ p l M _ l is the noise 
eigenvalue matrix. Furthermore, e p i spans the same signal subspace as that spanned by a(0 p ) . Thus, we 
have E^a(^) = 0 and SL H (0 p )E pn =0. The MUSIC estimates the DOA of the pth user from the 
highest peak of the following spectrum [3]: 

£*„(#) = max — — , p = \,2,---,P (Q\ 

where the largest maxima of S {0) is taken to be the estimate of the DOA of the pth user. 

As the angular spread A , of the Ath scattered ray from the pth source arriving at the array is 

EN 
h '-iYph a P h i n Equation (4) may be 

approximated as follows [14]: 

M (io) 

where P P = ^ J ^l l 7 P h^ P h /^^Tph an & & P = d& p /d0 p is the first derivative of the steering vector. 

Equation (10) is a linear model and regards as the generalized array manifold (GAM) model with 
scattering, which is the superposition of the pth source steering vectors and their derivatives. Then, the 
received data vector x(t) can be expressed as: 

x(0 = Us(0 + n(0 (11) 

Therefore, the pth user's despread and sampled array vector signal y p can be represented as: 

p 

y p =xc p =Ls p u(0 p )+ X s q &(6 q ) + n p (12) 

?=i, q*p 

The correlation matrix of Equation (12) in the observation interval is given by: 

K =E{y p r p } = L\u(0 p )u H (0 p ) + x np I M (13) 

The eigendecomposition of Equation (13) can be expressed as: 

M 

K^VjpJprn =V pl e pl e H pl +V p „Ap„K> (14) 

m=l 
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With the GAM model, the MUSIC searching function of Equation (9) can be expanded into a 
MUSIC with local scattering. The MUSIC searching function derived by the steering vector of the 
GAM model in Equation (10) , which is termed as GMUSIC and can be expressed as: 

(#,/?) = max — } _ - _ , p = l, 2, P fi^ 

" y ,H ' o.p \Z s H A H (0)E pn E H pn Mm { ' 

where A{6) = [a(#) d(0)] and % — [1 pf . The parameter of p can be estimated by using the 
method of [3]. And, the E pn in Equation (15) is the noise subspace, which is the eigendecomposition 
of the autocorrelation matrix derived by the linearized steering vector in Equation (10). We can 
estimate DOA by searching the array manifold for direction vector and the largest peak of the function 
GSp(8, p) denotes the DOA of the pth user. 

2.3. Problem Formulation 

Recently, a much more computationally efficient and more accurate parameter estimating approach 
via polynomial rooting method has been proposed to improve the spectral searching approaches for 
reducing the computational load. Due to the fact the scanning vector A{6) of the GMUSIC is not 
Vandermonde structure, GSp(6, p) cannot be implemented by using the polynomial rooting approach. 
The computational complexity of the classical subspace approach of MUSIC estimator is about 12M 
complex multiplications (CM) for computing EVD of a M x M dimension matrix. Let the search 
number be B. Therefore, the total required CM for computing Equation (15) of the GMUSIC estimator 
is about 12M 3 + BM 2 . The performance of the abovementioned spectral searching approach is 
governed by the scanning grid size and the number of search grids while implementing the 
high-resolution DOA estimation. A major problem of GMUSIC type algorithms is the heavy 
computational load incurred with spatial spectral search when a smaller grid size is adapted. However, 
it is time consuming and the optimum search grid size is not clear. Smaller grid size can improve 
estimate accuracy, but the required computational load also becomes relatively larger (e.g., if grid sizes 
are set to 1°, 0.1° and 0.01°, then the required search numbers B = 181, 1,801, and 18,001, 
respectively.). With the number of particle population P N and the maximum number of iterations & max , 
the computational complexity of our APSO approach is about P N x A: max . Hence, the computational 
complexity of APSO is smaller than the GMUSIC. Therefore, in order to reduce computational load, 
this paper investigates the feasibility of applying the PSO to replace the spectrum searching approach. 
As a result, the proposed PSO-based searching GMUSIC estimator does not increase the 
computational complexity, but significantly reduces the searching process requirement compared with 
the spectrum searching GMUSIC. 

3. The Proposed PSO Based Searching Algorithms 

Due to the fact that the performance of the abovementioned spectral searching GMUSIC estimator 
is governed by the scanning grid size and the number of search grids while implementing the 
high-resolution DOA estimation, it is time consuming and the search grid is not clear. In this section, 
we present PSO based searching approaches by maximizing the fitness function of GMUSIC at each 
iteration. 
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3.1. HPSO-GMUSIC Estimator 

In order to reduce the scanning accurate angle problems of high computation cost, we use the PSO 
to replace scan in the ULA. For time-space signal processing, each particle can be treated as a point in 
the P-dimensional problem space and a swarm consisting of D random particles and then searches for 
best position (solution) by updating generations until exceeding the limit of iteration number. The 
position of the z'th particle is represented as 8, = [0,i, 8,2, 8y>], the rate of the position change 
(velocity) for z'th particle is represented as v,- = [v,i, Vn, -, v,/>]. The best previous position of the z'th 
particle, which gives the best fitness value, is recorded as = [pn, pa, -, Pip]. And, the index of the 
best particle among all the particles in the population is represented by p g = [p g i, p g 2, Pgp] and called 
global best location. In every iteration, the local bests and global bests are determined through 
evaluating the fitness values of the current population of particles. Two factors characterize a particle 
status on the search space: its position and velocity. The P-dimensional position for the z'th particle in 
the Mi iteration can be denoted as G^/c) = [Qn(k), 6 i2 (/c), Q,?(k)]. Similarly, the velocity for the z'th 
particle in the Mi generation can be described as \i(k)= [v,i(/c), v,- 2 (/c) Vip(k)]. The velocity and the 
position of the z'th particle at the (k + l)th iteration for i = 1,2, D and k= 1,2, /c max are updated 
according to the following equations: 

yXk + ^-wi^y^ + c.r^^ip^-Qm + ^rA^^iV^-W) (16) 

e < (*+i) = e < (*)+v / (*+i) (17) 

where ® denotes element- wise product and the positive acceleration constants c\ and ci are two 
positive constants named learning factors. ri(/c) and r\(k) are P-dimensional vectors consisting of 
independent random numbers uniformly distributed between 0 and 1 , which are used to stochastically 
vary the relative pull of p, and p g in order to simulate the unpredictable component of natural swarm 
behavior. The inertia weight w(k) based on the linear decreasing strategy is considered critically for the 
convergence behavior of PSO [10], which is selected to decrease during the optimization process. Thus, 
PSO tends to have more global search ability at the beginning of the run while having more local 
search ability near the end of the optimization. Given a maximum value w max and a minimum value 
Wmin, w(k) is updated as: 

Mk + 1) = w max - Wmax " Wmm /V (18) 

max 

where /c max is the maximum number of iterations, w max and w m m are chosen as 0.9 and 0.4 respectively 
in this work as in [10]. Because the PSO particle represents a series of priorities that range from -90° 
to 90°, all parameters of the P-dimensional particle positions, either initialized or updated during 
search, must be confined within [6 m m, #max] = [~90°, 90°], avoiding infeasible particle positions that 
can lead to slow PSO searches. The new particle position is calculated using Equation (17). 

A particle position vector is converted to a candidate solution vector in the problem space through a 
suitable mapping. The score of the mapped vector evaluated by a GMUSIC function is regarded as the 
fitness of the corresponding particle. For time-space signal processing, the DOA GMUSIC estimation 
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problem in this paper is to find the z'th particle position or estimate angle 9 ip (k) of the pth user to 
maximize the following fitness function: 

fitness, (k) = „_„ ~* - 5 _ , » = 1, 2, P HQ^ 

At the end of iteration, global best location p gp (k max ) is the final estimated angle defined as 

@p = Pgp(^max)- The final global best position p g (/c max ) is taken as the GMUSIC estimates of users. 

As a result, the PSO with linear decreasing inertia weight, particle velocity limitation, and particle 
position clipping is termed as hard-constraint PSO GMUSIC (HPSO-GMUSIC). Unfortunately, the 
HPSO-GMUSIC is not so facile to implement for its overhead on computing the search space mapping, 
particle velocity limitation and particle position clipping. In this paper, we present a more feasible and 
efficient modified PSO algorithm for the HPSO-GMUSIC estimator. 

3.2. The Proposed Multiple Adaptive Inertia Weight 

The inertia weight is critical for the performance of PSO, which balances global exploration and 
local exploitation abilities of the swarm. In other words, the inertia weight w(k) is employed to control 
the impact of the previous history of velocities on the current velocity, thereby influencing the 
trade-off between global and local exploration abilities of the "flying points". A large inertia weight 
tends to facilitate searching new area and global exploration. Conversely, a small inertia weight 
facilitates local exploitation in the current search area. Suitable selection of the inertia weight provides 
a balance between global and local exploration abilities and thus requires less iteration on average to 
find the optimum. During the search every particle dynamically changes its position, so every particle 
locates in a complex environment and faces different situation. Therefore, each particle along every 
dimension may have different trade off between global and local search abilities. According to [15], it 
has been shown that the performance of the PSO algorithm with linearly decreasing inertia weight has 
the ability to quickly converge, but the PSO may lack global search ability at the end of a run and may 
fail to find the required optima in cases when the problem to be solved is too complicated and complex. 
But to some extent, this can be overcome by employing a self-adapting strategy for adjusting the 
inertia weight. In this subsection, inertia weight is dynamically adapted for every particle along every 
dimension. A measure, individual search ability, which characterizes the faced situation for every 
particle is defined. Basing on this measure, the particle could decide to whether to increase or decrease 
the values of inertia weight by means of the transfer function. The fine strategy of dynamically 
adjusting inertia weight could lead to improvement in performance of PSO. 

For velocity updating, the last two parts of Equation (16) can be viewed as the accelerations parts 
and can be defined as a ip (k) = c\-rand{*)-{p ip {k) - 9 ip (k)) and a gp (k) = crrand(*)-(p ip (k) - 9 ip (k)), where 
rand{*) is independent random number with uniformly distributed between [0,1]. So we could consider 
that the particle is moving with velocity of v ip (k) and acceleration of a ip (k) and a gp (k). But, a gp (k) is the 
dominant term for improving convergence rate. Suppose that the mass of the z'th particle in the pth 
dimension is normalized to 1 kg. According to the principle of mechanics, a gp (k) =f ip (k), where fi P {k) is 
an outside force, which is put on the particle comes from the "pulling" of pi P {k) and p gP {k). For ULA, 
the DOA is a one-dimensional searching problem. In order to make the particle fly towards optimal 
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region quickly, v ip (k) should turn to the direction off ip (k) as soon as possible. Define the error between 
Pi P (k) and 0 ip (k) is z ip (k) =p gp (k) - 0 ip (k). Let m ip (k) = \z ip (k)\/0 mSiX which is a number between 0 and 1. 
The velocity update problem of the z'th particle on the pth dimension can be divided into two classes: 

(1) Firstly, v ip (k) and z ip {k) are in the same direction, z ip (k)/v ip (k) > 0. If m ip {k) is relatively large, it 
means the particle is in the right direction, but the velocity is too small. Therefore, the particle 
needs to speed up, and the inertia weight w ip (k) needs to be set larger. If m ip (k) is relatively 
small, it means the particle has come to the location, that is near the optimal region. So the 
velocity of this particle should slow down and the neighborhood of the state should be searched 
carefully, and w ip {k) should be set smaller. 

(2) On the other hand, consider v ip (k) and z ip (k) on different direction, z ip (k)/v ip (k) < 0. If m ip {k) is 
relatively large, it means that the particle's state is far from the optimal region. So the particle 
needs to change its velocity as soon as possible, and the inertia weight w ip {k) needs to be set 
smaller. If m ip {k) is relatively small, it means that it's not urgent for the particle to change its 
direction on this dimension, and w ip (k) could be set a large value. 

Based on the aforementioned analysis, an adaptive inertia weight strategy is proposed and is shown 
in Table 1 . The individual normalized search ability of the z'th particle along the pth dimension is 
defined as m ip (k). It is noted that the value of w ip (k) is a function of m ip (k). To get a balance of global 
search and local search ability, w ip {k) cannot be too large or too small, thus w ip (k) is limited in the 
range of [w mm , w max ], which is like the process of normalization. Therefore, we used ju-law algorithm to 
achieve our strategy for every particle along every dimension, normalization of the m ip (k) to w ip (k) . 

The ju-law algorithm is a companding scheme used in telephone network [16]. Increasing the value of 
u, the dynamic range capability of ju-law can be improved and defined by: 

\og(\ + jUm.(k)) 
log(l + //) 

where m ip {k) and w ip (k) are normalized input and output values respectively, and the positive 
constant jU is the compression parameter. The reciprocal slop of this curve is given by the derivative 
of m ip (k) with respect to w ip (k) . For case 1, the value for small m ip {k) increases at the expense of the 
value for large nii P {k). To accommodate these conflicting requirements {i.e., a reasonable values for 
both small and large m ip {k)), a compromise is usually made in choosing the value of parameter ju for 
ju-law. The requirement of case 2 is opposite. Finally, based on our strategy and Equation (20), we can 
define multiple adaptive inertia weight w ip (k) as: 

mm , P K h 

z (*) (21) 



Hfc(*) = 



In Equation (21), w m in be added to avoid particles from stopping moving. The curves of w ip {k) with 
>v m in = 0 using different pi can be plotted in Figure 1. Then, we also investigate the sensitivity for 
APSO-GMUSIC with different values of fi However, it accords with our strategy for different //. Thus, 
//-law algorithm with \x = 100 is chosen. Note that for every particle in population, w ip (k) is unique and 
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can be computed individually. Therefore, the single inertia weight w ip (k) can be replaced by a multiple 
adaptive inertia weight w ip (k). The proposed APSO-GMUSIC seems to be robust to control parameters 
due to the intrinsic advantages of the algorithm and the separation of the problem-independent PSO 
kernel from newly introduced problem-specific features in our design for adaptive multiple inertia 
weight. Finally, the steps for implementing the APSO-GMUSIC are shown in Figure 2 and described 
in the list that follows. 



Table 1. Adaptive inertia weight strategy. 



Value of w ip (k) 


Value of m ip (k) 


Small 


Middle 


Large 


Directions of 

vtp (k) and z ip (k) 


Same 


small 


middle 


large 


Opposite 


large 


middle 


small 



Figure 1. Normalized output w ip (k) versus normalized input z ip (k)IO m?ai under w min =0. 
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Figure 2. Flowchart illustrating the main steps of the APSO-GMUSIC. 



*Initialization 

• Randomize particles' positions 

• Randomize particles' velocities 

• Search local and global best 
*Repeat for each iteration 

*Repeat for each dimension 
*Repeat for each particle 
*Particle Evaluation 

• According to fitness value according to ( 1 9) 

• Update local best position and global best position 
*Inertia Weight Updating 

• Compute w^A;) according to (20) and (21) 
*Velocity Updating 

• Compute v lp (k) according to (16) 
*Position Updating 

• Compute 9 ip (k) according to (17) 
*Next particle 

Test termination criterion 

• If meet, solution is final global best position g p 
*Next dimension 
*Next iteration 
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4. Computer Simulations 

Several computer simulations for illustration and comparison are presented in this section. We use 
a 6-element ULA with half wavelength for simulation. Consider an asynchronous CDMA system 
with PN codes of length 31 and BPSK modulation. The spatial signature u p is generated by 

N p = 30 independent identically distributed local scatters with uniform random angular spread of 
width 2A p . It is assumed that A p is uniformly distributed over the interval [~Ap, A p ]. Each scatter 

has equal power and is randomly, uniformly distributed over [0, 2n\ The value of parameter p is 
assumed to be known. Eight users are impinging on the array with random impinging angles of 
uniform distribution in the interval [-90°, 90°]. All signal powers are set to be 20 dBW. The additive 
background noise is assumed to be white Gaussian distribution with zero-mean and unit Watts of 
power. For GMUSIC, the search resolution (grid) is set to 1°, 0.1°, and 0.01° under the search 
range of [ _ 90°, 90°]. The root-mean-squared error (RMSE) is calculated in an average manner as 
RMSE = [(1/FP) Ylj=iYJp=i{Pip ~ fp) 2 ] 0,5 :. where F indicates the number of independent simulation 
runs, P is the number of users, 0j p is the pth estimate DOA achieved in the yth run, and 0 P is the true 
DOA of the pth user. The PSO parameters chosen for all the experiments are summarized in Table 2. 
All PSO algorithms start with random initializations and are terminated if the maximum iteration & max 
is reached or the global best particle position is not updated in 20 successive iterations. Every 
simulation result is presented after 200 data bits were processed and it is averaged by F = 10 3 
independent Monte Carlo runs with independent noise samples for each run. 



Table 2. Selected PSO Parameters. 









The number 
of particles 


The number 
of iterations 


Constraints 


Vfmax 




Inertia weight 
feature 


HPSO 


2 


2 


20 


20 


Hard 


0.9 


0.4 


Single and 
linearly 


APSO 


constraint 


1.1 


0.1 


Multiple and 
used jU-law 



Comparison results with other estimators, including the GMUSIC, HPSO-GMUSIC and 
APSO-GMUSIC with ju = 100 to DOA estimation error are presented. Figure 3 depicts the 
convergence in terms of DOA RMSE versus the number of iterations. As a result, the HPSO-GMUSIC 
requires more iterations to achieve convergence. Note that the proposed APSO-GMUSIC achieves fast 
convergence with the selected parameters, which means that it needs less iterations to approach the 
desired global extreme. Figure 4 shows the required number of calculating fitness function (B) versus 
number of particles. For the number of particles in the population, more particles may increase success 
in searching for optima due to sampling state space more thoroughly. However, more particles require 
more evaluation cost. The HPSO-GMUSIC needs more particles to approach the desired global 
extreme. It is confirmed that the proposed PSO-based searching approaches can reduce the 
computational complexity of the GMUSIC due to the searching process. As expected, this figure also 
provides a great improvement over the convergence rate on optimization problems. In fact, additional 
adaptive multiple inertia weight operation can improve the searching speed and RMSE performance 
further. Figure 5 presents the RMSE of DOA estimation versus varying angular spreads. We note that 
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the subspace-based techniques show serious degradation when faced with local scatters. Local 
scattering may be viewed as a form of model error and gives rise to the perturbation of the noise 
subspace. Again, these figures show that the proposed APSO-GMUSIC method yields significantly 
superior performance over the other methods in the presence of local scatters. For comparison, the 
result of GA-GMUSIC estimator is also provided. The same parameters of GA-GMUSIC estimator are 
used in [17]. Figure 6 shows the RMSE versus different SNR of the desired user for angular spreads 
2A P = 1°. For the low SNR case, all of methods may produce highly biased estimates. Clearly, with 
the compatible searching resolution, the APSO-GMUSIC can save the required number of searching 
grids and improve the RMSE performance, as compared with the other estimators. The GA-GMUSIC 
has a poor performance than HPSO-GMUSIC and APSO-GMUSIC. It is well known that premature 
convergence degrades the performance of GA and reduces the search ability [18]. In addition, it has 
been shown that the performance of the PSO algorithm with linearly decreasing inertia weight has the 
ability to quickly converge, the PSO may lack global search ability at the end of a run and may fail to 
find the required optima in cases when the problem to be solved is too complicated and complex [19]. 
But to some extent, this can be overcome by employing the proposed adaptive multiple strategy for 
adjusting the inertia weight. Finally, in Figure 7, we compare the RMSE performance against the 
number of active user P, given SNR = 20 dBW and angular spread 2 A p =1°. Basically, the RMSE is 
increased quite steadily with the increase of P. It can be observed that the APSO-GMUSIC obtain 
more performance improvement when the number of users P is reasonably increasing. Among them, 
the proposed APSO method achieves the lowest RMSE. 



Figure 3. DOA RMSE versus the number of iterations. 
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Figure 4. The required number of calculating fitness function (B) versus number of particles. 
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Figure 5. RMSE of DO A estimation versus angular spread. 
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Figure 6. RMSE of DOA estimation versus the input SNR of desired user for 2A P = V 
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Figure 7. RMSE of DOA estimation versus the number of users for SNR = 20 dBW and 2A P = V 
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5. Conclusions 

In this paper, we have presented a PSO based searching DOA estimation method, named APSO, 
which uses the GMUSIC searching function for CDMA signals. With the code-matched filter, the 
MAIs after code-decorrelation appear as noises for CDMA signals. The rooting MUSIC method is 
suboptimal in the presence of the noise and MAI and the GMUSIC cannot be implemented by using 
the polynomial rooting approach. However, the proposed techniques reduce the required search grids 
for the conventional spectral searching estimators. In the previous work on DOA estimation using PSO 
algorithm [10], hard constraints have been taken, during each iteration of PSO algorithm. In this paper, 
multiple inertia weight has been incorporated in PSO. In conjunction with a modified PSO for angle 
searching, the proposed approach can reduce the required computational load for the conventional 
spectral searching MUSIC estimator with local scattering. Moreover, the convergence of the proposed 
approach is much faster. Computer simulations have demonstrated the effectiveness of the proposed 
approach. Furthermore, a common drawback with MUSIC-like technique is that it is a suboptimal 
estimator and tends to suffer from low performance due to low SNR, small sample size, and correlated 
sources. Thus, how to develop a cheaper way and work well with correlated or even coherent sources 
will be important in the future work. 
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